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We propose a method for the numerical computation of microcanoni- 
cal expectation values — i.e. averages over energy eigenstates with the same 
eigenvalue — without any prior knowledge about the spectrum of the Hamil- 
tonian. This is accomplished by first defining a new Gaussian ensemble to 
approximate the microcanonical one. We then develop a "lattice theory" for 
this Gaussian ensemble and propose a Monte Carlo integration of the expec- 



ts 

■ tation values of the lattice theory. 
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I. INTRODUCTION 

The microcanonical ensemble is considered to be the fundamental ensemble of statistical 
physics. For example, the use of the canonical Gibbs ensemble is usually justified by showing 
that its expectation values coincide with the microcanonical ones in the infinite- volume limit. 
However, since an application of the microcanonical ensemble requires detailed knowledge 
about the energy levels of the system, it is seldom possible to use it in practice. The canonical 
ensemble, on the other hand, has a well-behaved path-integral expression easily extended 
OS . to the gauge field theories — therefore it has become the standard ensemble of quantum 
^ ! statistics. 

Nevertheless, there are situations where the easiest alternative does not work. For in- 
stance, direct applications of the grand canonical ensemble have not been able to reproduce 
all the results of the relativistic ion collision experiments!] There are two possible explana- 
tions for this shortcoming of the canonical ensemble: either the particle gas created in the 
collision does not reach thermal equilibrium before exploding into the final state hadrons, 
or the system is too small to be handled by the canonical methods. 

Recent calculations! using the grand canonical ensemble with finite-volume corrections 
have, however, succeeded in describing most of the particle abundances in heavy ion colli- 
sions. This suggests that at least the final state hadron gas will thermalize, but it also shows 
that the finite-volume effects are prominent in these systems. 

These results point into the direction that the microcanonical ensemble — which is optimal 
for describing isolated ergodic systems with small quantum numbers — should be used for 
getting quantitative information about the properties of the quark-gluon plasma possibly 
created in the relativistic hadron collisions. The continuum path integral formulation of 
the microcanonical ensemble and the corresponding perturbation theory have already been 
developed both for scalaiH and gauge field theories^ but since the phenomena involving 
hadrons happen in the low energy, strong coupling regime, non-perturbative methods are 
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also needed. However, no general numerical treatment is presently available to compute 
microcanonical expectation values without a priori knowledge about the energy spectrum 
even in the quantum mechanical context, let alone in a quantum field theory. In this 
paper, we aim at a development of such a method in the spirit of lattice field theories 
for a class of quantum mechanical Hamiltonians. The generalization of these results to the 
more interesting case of gauge field theories, such as quantum chromodynamics, is a subject 
of further work. 

The paper is divided into three parts: the first section introduces a new statistical ensem- 
ble, called a Gaussian ensemble, and establishes its close connection to the microcanonical 
one. In the second section we will derive lattice integral kernels for a numerical evaluation 
of the Gaussian expectation values. Due to the complicated nature of these kernels, it is not 
immediately evident how Monte Carlo methods can be used in this computation — this will 
be the subject of the rest of the paper. 

II. THE GAUSSIAN ENSEMBLe! 

Consider the operator defined by 



V2ne 



1 I H — E 



(1) 



where H is the Hamiltonian. This operator is clearly bounded, self-adjoint and positive. 
Let us assume that the spectrum of the Hamiltonian is discrete with eigenvalues E n , each of 
which has a finite number of eigenvectors Q n ,k, k — 1, . . . , K n , where K n is the multiplicity 
of the eigenvalue E n . With rather loose assumptions on the eigenvalues (lim^oo n s /E n = 
for some s > will suffice if the multiplicities n n diverge at most polynomially in n) p e is 
trace-class and therefore defines a sensible statistical ensemble via 

(A) = — ^— ^ — , for any observable A. (2) 

Note that if the Gibbs ensemble is well-defined for temperature 2e 2 , i.e. Tr exp(— -^H) < oo, 
then E n — E < 1 only for finitely many n, k and, since for the rest (E n — E) 2 > E n — E, f> e {E) 
is then trace-class. This means that the Gaussian ensemble is well-defined for all systems 
for which we can use the canonical one. 

Intuitively, the above density operator describes a system that was initially prepared into 
a Gaussian energy distribution E±e, but which has then stabilized under ergodic evolution 
in isolation. This, with the formal relation lim £ ^o Pe(-E') = 5(H — E), will give the physical 
motivation for using this ensemble to approximate the microcanonical one. However, it 
can also be argued that the Gaussian ensemble is even better suited for describing real 
experimental situations than the microcanonical ensemble, since it naturally incorporates 
the energy resolution of the experiment through the parameter e — in any real experiment 
the total energy of the system (e.g. of a particle) cannot be known exactly, but only within 
a certain accuracy. We will now turn to the mathematical justification of these formal 
statements. 
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Using the earlier notations, we can express the trace in (|2|) as 



Tr(pU) 



V2 



!7T£ 2 ^ 

n,k 



cxp 



{^n,k\A\fl njk ) . 



Since 



lim exp 



0, when E ^ E n , 



the partition function clearly satisfies 

limTrpe = Y] n n S(E - E n ) = Z mi 

e— >0 z — ' 



microcan- 



In the evaluation of the limit of the expectation values, we need the result 

1 f 0, if \E„i — E\ < \E„ — E\ for some n' 



E ex P 

n'k' 



(E n ,-E)*-(E n -E)* 
2p 



(En' K n'$\E n ,-E\,\E n -E\) \ otherwise 



(3) 



(4) 



(5) 



(6) 



In other words, this limit is zero if E n is not the eigenvalue nearest to E, it is (k„ + k„') _1 if 
both E n and E n i are nearest eigenvalues (i.e. if E lies exactly in the middle of the segment 
joining E n and E n > and no other eigenvalues are on this segment) and it is k" 1 if E n is a 
unique nearest eigenvalue. Let us use notation M n {E) for the sum En' K ri&\E n ,-E\ ) \E n -E\- 
Using the formula (H) for the traces and (0) to compute the limit, we arrive at 



exp (E n -E) 



n,k 



Tr f) e \Jl-ne 1 



(^ra,fc|^4|^n,fc) 



M no (E) 



^2S\ En -E\,\E no -E\(^n,k\ A \^n,k), 
n,k (7) 



where uq is the index of an eigenvalue nearest to E. Since M no (E) is the number of terms 
in the last sum, the last expression is nothing but the average of the expectation values of 
A over the energy eigenstates nearest to E. 

Therefore, whenever E coincides with a point in the spectrum, the limit e — > will give 
the microcanonical expectation value. On the other hand, if E does not belong to the energy 
spectrum (in which case the microcanonical ensemble is in principle ill-defined), then the 
microcanonical result corresponding to the nearest eigenvalue is obtained. The only values 
of E giving non-microcanonical limits are those lying exactly in the middle between two 
eigenvalues, but even then the result is an expectation value of a uniform distribution over 
two energy eigenvalues. 

In the above, we tacitly assumed that it is possible to exchange the order of the limit and 
the sum over n. This operation is in fact allowed if the sum in (|]) converges absolutely for 
any e, which for positive A is tantamount to the assumption (A) e < oo. We will now prove 
this by showing that if (|3|) converges absolutely for e > 0, then the series in (|7|) converges 
uniformly for < e < Eq. 

Let n be, as before, the index of one of the eigenvalues nearest to E. Since exp( — is 
an increasing function of e for C, e > 0, we immediately see that 
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E ex P 

n'k' 



( E n' -E) 2 -(E n -E) 2 



2e 2 



< exp 



{En - Ef - (E no - E) 
If 2 



(8) 
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for all n such that \E n — E\ > \E no — E\. In fact, this inequality extends also to 
the case \E n — E\ = \E no — E\, since the ratio on the left hand side is always less 
than one. Therefore, the absolute value of each of the terms of the series in (^) is less 
than exp[(E no — E) 2 /(2el)]exp[—(E n — E) 2 /(2el)]\(Q ntk \A\Q ntk }\, which again form an e- 
independent sum that is convergent by assumption. Thus the series in (^]) satisfies the 
Weierstrass M-test and is uniformly convergent. 

The uniform convergence of V2ne 2 Ti p e is similarly established by the Weierstrass M-test 
using the monotonicity of exp[— (E n — E) 2 /{2e 2 )}. 



III. DERIVATION OF THE LATTICE ENSEMBLE 

In this section we will derive a "lattice approximation" of the Gaussian ensemble for 
quantum mechanical Hamiltonians of the type H = ^—+V(x), where V is a suitable potential 
function. For example, if V were a (non-constant) polynomial bounded from below, then H 
would be essentially self-adjoint and the operator e~^ H would be trace-class with the usual 
Feynman-Kac path-integral kernel. Also, for any polynomial and for all bounded functions 
A, the canonical expectation values would then be well-defined and could be computed by 
using the formulaBEl 



Tr A{x)e-f> B ] = Jim / ^/^ ) exp ij>(**-i - x k ) - £ £ ^ + V(x k ) 

J ^ > L k=i fe^i ^ 



(9) 



where x = and /3 > 0. 

In the following we will need the analytic continuation of this result to the region Re j3 > 
0. If V is a second-order polynomial, i.e. if we consider a harmonic oscillator, then this 
continuation can be computed exactly and it is seen to be given by @. There is no apparent 
reason why the same continuation should not work for any polynomial, but we have not yet 
been able to construct a proof of this generality.^ We will therefore simply assume that V is 
a potential such that (§) holds for all Re (3 > 0. 

It will be useful to define a few auxilary quantities for the computation of the lattice 
expectation values. Firstly, the relevant energy parameter of the lattice Gaussian distribu- 
tion will turn out to be E/e and we will use a separate notation W for this dimensionless 
quantity. Secondly, we will need the following three "action-like" quantities: 



1 N / - \ 2 i N 

and, finally, S(p,x;e) = jj^J2k=i ~^ A comparision with (§) reveals that S is 

nothing but the usual finite-temperature lattice action for temperature T = e. 

In order to keep track of all the mathematical details in the derivation of the lattice 
ensemble, we need to use a temporary regulator to provide for absolute convergence of the 
integrals. In the following we have used exp(— (3H) for this purpose — the potentials have 
been chosen so that this is legitimate. 
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Since the energy spectrum is bounded from below, we can use this lower bound to give 
a majorant in the M-test and hence to justify the formula 



Tr (p £ A) = lim , 



: exp 



2e 2 



{E n - Ef - (3E n 



(fin,k | A\ ^n,k) ; 



for A = A(x) as before. The Fourier transform of a Gaussian distribution will now give the 
integral representation 



da 
2tt' 



-\e 2 a 2 +iaE-iaR 



V2 



(12) 



which, by using the bound Tr 
rem, will lead to the result 



\A\e P H < oo in the Lebesgue dominated convergence theo- 

00 da 



Tr (f> £ A) = lim / — e 



(3^0+ 



-±£ 2 a 2 +iaE 



2tt 



Tr 



-(p+\ a )H A 



(13) 



The trace in the integrand of ( |13"D has now a path-integral expression given by (^). In- 
serting the path-integral formula, referring once more to the dominated convergence theorem 
to move the "continuum limit" over the a-integration and applying the Fubini theorem and 
eq. (|12"D will then yield 



Tr (p E A) = lim lim 

VF ' 13 N 



d N xd N p A(x ) 
(2tt) n 



exp 



i^Pkixk-x - x k ) - (3eS - \{S - Wf 



If we now use the equality 



PeS - - {S - wy 



~(S -W + f3ef + ^(3 2 e 2 - f3eW 



(14) 



(15) 



and take the limit of the last two iV-independent terms, the effect of using the regulator will 
reduce to a shifting of the energy E by —j3e 2 . But since the series in (|TT1) clearly converges 
uniformly in E and thus gives a continuous function in it, the limit (3 —>■ + will amount to 
setting 13 = in (0). 

Therefore, we have proved that 



Tr(p e A)=lim / d N x A(x )K N (x; E, e), 

N—*oo 



(16) 



where the kernel is given by 

d N p giEPkOzfc-i-Zfc) 



K N (x;E,e) 



(27T 



V2 



exp 



2mNe 



Si(x;e) - W 



(17) 



Most of this integral can be evaluated by choosing spherical coordinates with the iV-axis 
pointing in the direction of the vector (xk-i — x^- Let us denote by y the rotated iV- 
coordinate and let us gather the remaining components into a N—l- dimensional vector q. 
Then the kernel can be rewritten as 



K 



N 



d N ~ l q 
(2tt) n - 1 



.2tt 



V2 



TVS'' 



exp 



1 / y 2 + q 2 

2 \2mNe 



+ Sx - W 



(18) 
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Next we will choose polar coordinates in the g-space and rescale both q and y by the factor 
1 / \ / 2mNe. Since the integrand is independent of the angles and the area of the unit sphere 
in IR n is given by 2vr n / 2 r(n/2)- 1 , we will then arrive at the formula 



N 



2(mNe)~ f°°, N _ 2 x 







i2y^S~ -±(y 2 + q 2 + S 1 -Wy 



J(19) 



where we have made use of the definitions in (|T0|). 

It is possible to re-express this Fourier transform as a one-dimensional real integral. If 
we choose polar coordinates in the (y, q) -plane, utilize the following integral representation 
for the Bessel functions of the first kind,! 

jyesm N ~ 2 9e ivcose = r (^7T^) V^{2/v)^ J^v), (20) 



and, finally, denote r = y 2 + q 2 and v = ^j^, we have 



wy 



(21) 



The integral representation fl2T|) can now be used to derive a series expansion for the 
kernel. Since 

the expansion of the exponential exp[— r{S\ — W)] and a subsequent term by term integration 
will give 

K _ (mNe\% exp [-1 ( Sl - W) 2 ] " " (-S )* (W - ^)"" fc 2^r(^) 

N { 2n ) V2^ k] n^+i+k) ' (23) 

which, however, has no immediate applications for evaluation of the kernel, because the 
convergence of the sum will be too slow for large values of the parameters. 

So far, the best way to evaluate the kernel for large enough range of parameter 
values has been a numerical integration of eq. (ET|). This is mainly due to the well-developed 
methods of computation of the values of the Bessel functions. Similar methods — recursion 
relations, integral representations, etc. — will likely offer more efficient ways of evaluating 
the kernel as well. However, since the numerical integration was sufficiently fast for our 
purposes, we do not pursue this point further here. 



IV. NUMERICAL EVALUATION OF THE LATTICE EXPECTATION VALUES 

The purpose of this section is to enable a Monte Carlo integration of the lattice expec- 
tation values in ([16]). From the experience in the lattice quantum field simulations, this 
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should be possible provided that the following condition is met: the Monte Carlo is used 
for a generation of a probability distribution exp(— S"[x]), where S[x] is a local action. The 
locality, which means that the action is a sum over a local "action density" , makes it possible 
to use local updating algorithms and in this way to achieve a fast enough generation of the 
required probability distribution. In the following we will derive an exponential bound for 
the kernel, but with an "action", which is not local. However, as explained at the end of 
this section, a certain linear approximation for the exponent will be a local quantity, using 
which it is still possible to benefit from the local updating algorithms. 

Letting r = q 2 in (|Tj|) will give the starting point for the approximations: 



K 



N 



(mNe) 



(2.)fr(^)V2^¥ 



oo poo 

drr 2 



dyexp 



i2yy%--(y 2 + r + S 1 -Wy 



J(24) 



Next we need to rely on the lemmas proved in the Appendix. First on Lemma |A.1| , which 
states that for all real A, the y in the integrand of ( p4|) can be replaced by y+iX, i.e. that 
the double integral is always equal to 



OO /"OO 
N=3 

drr 2 



dyexp 



i2y^/S~ -2X^/S~ -^ (y 2 -X 2 + r + S 1 -W + i2Xy) : 



(25) 



The absolute value of this quantity is less than 



dy exp — 2X>y S + 2X y / drr 2 exp 

'0 



1 -(r + y 2 -X 2 + S 1 -Wy 



(26) 



Now we are ready to apply Lemma |A.2| to get an upper bound for the second integral, 



dyexp -2Xy/ S + 2X 2 y 2 



2^r(M) 



f-exp[-b N ( y 2 -X 2 + S 1 -W)}, 



(27) 



with b N = v / 2r(^±i) /r(^p-) (for large values of N, bN can be approximated by x 1 
\py\ Since 



dy exp[-(b N -2X 2 )y 2 } 



71 



>N 



2A 2 



(28) 



whenever 2 A 2 < b n — L we can infer from (|27|) that for these values of A 



K N \ < 



AT 

fmNe\ 2 1 



exp 



-2A V^o + b N X 2 - b N {Si - W) 



(29) 



As a function of A, the right hand side has only one minimum which lies at A = v^So/frjv- If 
So < b 2 N {bN — l)/2, we can use this value of A to give the upper bound, and the exponential 
factor will become exp[— S /b N — b^(Si — W)\. 
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Let us now use notation Mjv for the critical value of So — for large N, Mn is approximately 

equal to If S > M N} the best allowed value is then given by A = \J i>N 2 1 = \/M N /b N , 

which yields the exponential factor exp [— {2^MnS$ — Mjv)/6w — b^{Si — W) 
To summarize, we have now proved that 



K N \ < 



fmNeX 



71 



exp 



Wb 



N 



So 



+ b N S 1 



x exp 



N 



(30) 



where 9 denotes the step function. As a curiosity, it should be noted that the term 
So/bN + bftSi corresponds to a certain "renormalized" finite temperature lattice action— 
renormalized in the sense that the temperature of the action depends on the size of the 
lattice: T{N) = e/b N . 

We propose that the lattice ensemble expectation values could be evaluated by using a 
Monte Carlo sampling to generate a distribution exp(— S^ N \x; E,e)) and thus allowing the 
evaluation of the ratiosH 



/ d N x A{x )K N {x; E, e) exp{-S^ N \x; E, e)) 
Jd N x K N (x; E, e) exp(—SW(x; E, e)) 



(31) 



where is obtained from the scaling of the kernel by the inverse of the right hand side of 
(30) and the "lattice action" S {N) is given by the formula 



S(N) = s -e(s -M N ) (VgW^y 



+ b N (S 1 - E/e) . 



(32) 



For local updating algorithms, the value of the action could be computed exactly at the 
beginning of each sweep and the linear approximation 



AS {N) 



-ASo + b N AS 1 



(33) 



could be used during the sweep. In the linear approximation, c = 1, if S < M N at the 
beginning of the sweep, and c = a/ M n / S otherwise. 

To check that the approximations used in the derivation of the upper bound are reason- 
able, we have shown a plot of values of the normalized kernel with N = 22 and W = 7 
in Fig. 1. Two general features are evident from the figure: firstly, there is a "gaussian wall" 
around the line Si = W, which means that the importance sampling should extend at least 
that far in the Si-direction. Secondly, although the scaling can be seen to be quite optimal 
for the region Sq < Mat, the scaled values drop rapidly after that point. The oscillations 
of the kernel, however, also begin at about the same values of <So and an oversampling of 
the oscillatory region is in fact only desirable — therefore, the use of this particular scaling 
function is qualitatively justified. 
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V. CONCLUSIONS 



We have defined a Gaussian quantum statistical ensemble for the approximation of micro- 
canonical expectation values of any system which is amenable for the usual Gibbs ensemble 
treatment. We have also derived a method for the computation of these expectation values 
by a Monte Carlo integration on a time-lattice, provided that the Hamiltonian allows an 
analytic lattice regularization. The generalization of these results for quantum field theo- 
ries and tests for the numerical efficiency of the lattice approximation are currently under 
investigation. 
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APPENDIX A: 

In this Appendix we shall prove the two lemmas that where used in Section [TV] to derive 
the exponential bound for the integral kernel. 

Lemma A.l Let f be defined by the formula f(y) = exp[iby — \{y 2 + C) 2 ] , where C and b 
are complex parameters. The integrals of f over lines parallel to the real axis are all equal: 

/DC POO 
dyf(y + i\)= / dy f(y), for all A G R. 
-oo J — oo 



Proof. Consider the integral of / over a countour following the perimeter of a square with 
corners at points ±R, ±R + iA. Since / is a composition of two entire functions, it is entire 
and the integral over this countour is therefore zero by the Cauchy theorem. Writing this 
integral in a parameter form gives 

[ R dyf(y)= [ R dyf(y + i\)-i\ [ dt(f(R + i\t)-f(-R + i\t)). 

J-R J-R JO 



But since for all < t < 1, 

\f(±R + iXt)\ < exp 



|6 J R| + |6A|-i( J R 2 -A 2 -|C|) 2 + i(2|A| + |C|) 2 



and the upper bound goes to zero when R goes to infinity, the integral over t vanishes when 
we take the limit R — > oo and the lemma follows. □ 
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Lemma A. 2 For all r > and CgI 

P CO 

/ drr T e-^ r+c f < a T e- b ^ c , 
Jo 

where a T and b T are positive constants given by the formulae 

VHXr + i) >/zr(fg) 
2 ^r(^)' ' T r(*i) ' 

Proof. Let us use the notation I T {C) for the integral J °°dr r T e~^ r+c ^ 2 '. The integral 
converges for r > — 1, when J T defines a function from the real numbers to (0, oo). We claim 
that the following three properties are satisfied for all C, if r > 0: 

(a) I t+1 (C) = tI t _ 1 (C)-CI t (C), 

(b) J T is differentiate at C and = -I T+ i(C) - CI T (C) = -rI T ^(C), 

Ir+l(C) > I AC) 



i T (c) " i T -i(c)- 

The lemma is then a fairly straightforward consequence of these results. 

Fix r > and define g(C) = ln/ T (C). Applying (|5|) twice shows that the first and the 
second derivatives of g exist and that they are g' = —rI T -\jI T and g" = —r^^j^-^f 1 — 1). 
But from @ we now deduce that g" < and therefore g is a concave differentiable function. 
This means that every tangent line to g is an upper bound for the function: for every fixed 
DgK, g{C) < (C — D)g'(D) + g(D). Exponentiation of both sides of this inequality gives 
now 



I T (C) < / T (D)exp 



- T ~[jDT ic - D \ 



which is of the desired form if we just define b T = tI t _i(D) / I T (D) and a T = I T {D) exp[b T D}. 

It now remains to show that the explicit formulas given in the Lemma for a T and b T 
are valid for some D. Choose D — 0, when a change of variables from r to \r 2 shows 
that / r (0) = 2 i 2 jL r( 1 ji). On the other hand, an application of (|J) reveals that for D — 0, 
b T = J r+1 (0)// T (0), which then simplifies into the desired formula. Since now a T = I T (0), 
the form for a T is a consequence of the Gamma-function "doubling formula": r(^ii) = 

2^(1) r(z). 

Therefore, it is necessary only to show that the propositions (a)-(c) hold. (|a|) is obtained 
by a simple partial integration, since r exp[—|(r+C) 2 ] = — 4- exp[— |(r+C) 2 ]— Cexp[— |(r+ 
C) 2 ] and the boundary terms vanish for r > 0. 

Consider next the finite difference term, I T (C + h) — I T {C). This is equal to 
J °°rfr r T e~2( r+c ') 2 {exp[— \h 2 — h{r + C)\ — 1}, where the term in the braces is bounded 
by exp[|/i|r + \hC\] + 1. Lebesgue dominated convergence theorem now shows that it is 
possible to differentiate under the integral sign, which will prove the differentiability and 
the first equality in ([Bp. The second equality follows then by an application of (^). 

To prove (^), it is enough to show that J T+1 (C)/ r _i(C) — I 2 {C) > 0, for all C and r > 0. 
Expressing the products as two-dimensional integrals will yield 
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poo poo 

/ T+ i(C)/ T -i(C) - I 2 T (C) = dx dyx T y T - l (x-y)exp 

Jo Jo 



- l -{{x + Cf + {y + Cf) 



In the change of variables (x, y) — > (t,u) — (x + y,^r), the Jacobian equals \t and the 
integration region changes into (0, oo) x (—1, 1). This way the integral becomes 



poo pi 

J dt J du2- 2T t 2T+1 (u + u 2 )(l-u 2 y- 1 exp 



-U 2 (l + u 2 )-Ct-C 2 
4 v ; 



As the region of integration over u is symmetric, the antisymmetric part vanishes, and we 
arrive at the result 



I T+1 (C)I T ^(C)-I 2 (C) - 
which is always strictly positive. □ 



(// / du2~ 2T t 2T+1 u 2 (l-u 2 ) T - 1 exp 



--t 2 (l+u 2 )-Ct-C 2 
4 
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FIG. 1. Values of the normalized kernel Kn as a function of the variables So and Si from two 
different viewpoints. Here W = E/e = 7 and N = 22, which corresponds to v = 10, 6jv ~ 3.24 
and Mjv ~ 11.8. The maximum value at the peak near So = Mn, Si = W is approximately 0.67. 



